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We study cosmological expansion in F(R) gravity using the trace of the field equations. High 
frequency oscillations in the Ricci scalar, whose amplitude increases as one evolves backward in 
time, have been predicted in recent works. We show that the approximations used to derive this 
result very quickly breakdown in any realistic model due to the non-linear nature of the underlying 
problem. Using a combination of numerical and semi-analytic techniques, we study a range of 
models which are otherwise devoid of known pathologies. We find that high frequency asymmetric 
oscillations and a singularity at finite time appear to be present for a wide range of initial conditions. 
We show that this singularity can be avoided with a certain range of initial conditions, which we 
^vq , find by evolving the models forwards in time. In addition we show that the oscillations in the Ricci 

scalar are highly suppressed in the Hubble parameter and scale factor. 
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I. INTRODUCTION 



It is now generally accepted that the Universe is currently undergoing a period of accelerated expansion. The 
most popular approach to modeling the current epoch is to postulate a new energy component of the Universe, 
dark energy, which has negative pressure and hence drives the acceleration. Possible dark energy candidates include 
the cosmological constant which has well known fine tuning issues, quintessence fields (see for example Q and 
I ■ references therein) and elastic dark energy 

An alternative approach is to modify the theory of gravity, and there have been a number of attempts to construct 
t-^ . late time accelerating solutions to the gravitational field equations by modifying gravity at large distances/late times. 
One such class of models are found by considering a gravitational action that contains an arbitrary function F(R) of 
the Ricci scalar. The first modified gravity model to be considered introduced quadratic terms into the gravitational 
action 0, IE IE 0| ■ Such models can give rise to early-time inflationary solutions to the gravitational field equations, 
where the Hubble parameter is initially H ~ M p i, and slowly rolls towards a stable Minkowski vacuum state. 

Modified gravity models which yield late-time acceleration have been considered recently. An important example 
is the CDDTT model Q, for which F(R) = R — jjft/R. This function has a de-Sitter vacuum solution to the field 
f^i . equations, which might be associated with the current epoch of the Universe with a suitable tuning of the mass scale 
p.. However, it has been shown that this model cannot satisfy local gravity constraints [9(, and also has an instability 
[Io| . In addition to the CDDTT model, a number of modifie d g ravity functions have been proposed which give rise 
to late time acceleration [lj], EE EE EE EE EE EE EE EE HE HE HE H3i an d some progress has been made in 
understanding general properties of these models, for example the behaviour of metric perturbations [2J, |25[ and the 
cosmological evolution of generic F(R) functions [IE HE HI] have been considered. 

Any F(R) model must satisfy certain conditions in order to exclude the possibility of ghost degrees of freedom and 
other instabilities. These conditions are F'(R) > 0, which ensures that there are no ghost degrees of freedom in the 
model, and F"(R) > 0, which prevents instabilities from arising in the early Universe [29|]. It seems also sensible to 
impose F(Q) = 0, which ensures that there is no cosmological constant, and F(R) — » R for large R, since we want any 
modification to gravity to only become significant at late times. In addition to these conditions, it has been shown 
that if a standard matter era (for which a oc t 2 ^ 3 ) is to be a fixed point then further conditions must be satisfied [2(|. 

Recently, a number of models have been studied in the literature which satisfy the above constraints, at least for 
R > 0. Specifically the following functions have been considered [29l. [30l. [3l| . 
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Fnss(R) = R-^ v 7\ 2n , (1) 



F AB (i?) = | + ^log 



, . 2bR\ , , . , ( 2bR 
cosh — tanh b smh 



(2) 



where i? V ac is the current vacuum curvature of the Universe, b, c are dimensionless constants, and n > 0. In ref . [29j| . 
the n is defined differently to ([T|), and in ref.[3l[ the precise formula was slightly different to (JTJ) . However, both 
models considered in (29l. [3l| have the same expansion in the limit R 3> -R va c, which will be the important point in 
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this paper. The model -Fhss represents a modification to General Relativity that is a power law in R, whilst Fab 
contains exponential corrections. These models both have late-time accelerating epochs, and can satisfy local tests of 
gravity. They can also satisfy the conditions discussed in ref . (20l| . 

Although the models (|1I2[) possess many desirable features, an additional issue has arisen recently with regards to 
their suitability as viable gravitational models [3l|, • Specifically, it has been suggested that if we look for solutions 
to the gravitational field equations that are perturbations around known General Relativistic solutions, then the Ricci 
scalar oscillates over very short timescales, and the amplitude of these oscillations will increase without bound to the 
past [3l|. If true, this presents us with a number of issues, amongst which are that the small perturbation will become 
much larger than the General Relativistic solution at some time in the past, violating the perturbative assumption, 
and R will become negative at some finite time. Such phenomena would most likely lead to significant problems in 
reproducing the standard Cosmology. 

The aim of this paper is to study the two model s (I1I2D in detail. We find that since R oscillates over such small 
timescales, the linearized analysis presented in refs.[3l|, l33| will not necessarily yield a solution that is indicative of 
the full solution for all times. We then solve the gravitational field equations numerically and semi-analytically and 
show that the Ricci scalar undergoes asymmetric oscillations, and will drift away from the General Relativistic limit. 
We note that non- linear oscillations have been considered in F(R) models previously, see ref.[34|. 

II. MODIFIED GRAVITY FORMALISM 

In this section we consider the modified gravity action 

S = J V^9d 4 x (M^F(R) + £ m ) , (3) 

which yields the following field equations 



1 T 



F'{R)R^ - ^uF(R) + [g^n - V P V„] F'(R) = (4) 



RF'(R)-2F(R) + 3nF'(R) = ^-, (5) 

where primes denote differentiation with respect to R, is the energy momentum tensor and is the reduced 
Planck mass. 

The HSS and AB models can be expanded as 



F(R)^R-^+ x (R), (6) 

for R > i? V ac, where i? V ac acts as a small 'cosmological constant' (although there is no true constant in these models 
since globally F(0) = 0). x(R), x'(R) and x"(R) are all small functions of R, in the sense that the dimensionless 
parameters x(R)/R, x'(R) and Rx"{R) all satisfy x(R)/R < 1, x'(R) < 1 and Rx"(R) <C 1 for i? > i? vac . For the 
HSS and AB models, we have 



2n+l 

Xhss = -jpr> Xab = £ab exp (-R/cab) , (7) 

where i? V ac is the current vacuum curvature of the Universe, and cab = -Rvac/46 and enss = -Rvac/(2c) 1 ^ 2n+1 ^ are 
parameters that are smaller than i? V ac- These expansions are valid for R ^> enss and R ^> eAB for the HSS and AB 
models respectively. For the rest of this paper, we will drop the subscripts ab and hsSj and use e = enss or eAB and 
X — Xhss or Xab, which should be obvious in context. We will explicitly state which model is being studied in each 
section. Finally, we note that we will frequently make use of the notation Rqr, which is the General Relativistic 
solution to the field equations Rqr — ~T /M^ v In particular during the matter era we have Rq-r = 4:/3t 2 , and during 
the radiation era Rgr oc i~ 3 / 2 . 
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III. 



PERTURB ATIVE ANALYSIS 



In this section, we review the perturbative analysis that has been used in refs. [3lL 133 . [33j to derive approximate 
solutions to the gravitational field equations for the AB and HSS models. The full field equations are a set of non- 
linear, fourth order differential equations for the scale factor, and solving these equations directly is difficult (although 
see for example ref. 35[ for an attempt to do so for the CDDTT model). The approach taken in this section is to 
look for solutions to the field equations that are perturbations around known General Relativistic solutions, that is 
we look for a solution to ([5|) of the form R = R G r + SR, where SR <C Rgr is some small perturbation. 

This approach involves linearizing the equation in SR. In doing so, the following expansions are used 



F(R) 

F'(R) 
F"(R) 



R GR + SR- 



R v 



+ x(Rgr) + x'(Rgr)SR + O(SR) 2 



In which case §5§ becomes 



1 + x'(Rgr) + x"(Rgr)SR + 0{SR) 2 , 
x"{Rgr) + x"'(Rgr.)SR + O(SR) 2 . 



3 X "DSR + 6 X '"V a RGRV a SR -5R = a(R GR ) 



(8) 
(9) 
(10) 



(11) 



to linear order in 8R. Unless otherwise stated, in (lll[) and for the remainder of this paper x is taken to be a function 
of i? G R only, x = x{Rgk)- The function a(R GR ) is 



v(Rgr) = 2x - Rgrx' ~ 3x'"(Vi? Gfl ) 2 - 3 X "aR GR , 



(12) 



which is a function of Rqr only and is small. The field equation (| 1 1 1) is an inhomogeneous, second order differential 
equation for SR. The solution is the linear sum of a particular solution to (|TT|) and the solution to the homogeneous 
equation 



3 X "DSR + Q x '"V a RGRV a SR -SR = 0, 



which can be written approximately as 



where a(t) is the scale factor. 



A. HSS model 



(13) 



(14) 



In ref.[31|, a solution to p4|) was derived for the HSS model. Using x = Xhss and e = chsS) the equation for SR 
reads 



2n+2 



SR + 3ff- 4 ^ 1)i?GR W R ™ 



Qn{2n + 1) 



1)e 2n+i 6R ~ °> 



(15) 



A WKB solution can obtained for the matter era, 



SR mat = i~ 3 "- 4 [At sin (A 2 t- 2n - 1 ) + A 3 cos (A 2 r 2n ~ 1 )] , 



(16) 



and for the radiation era 



«?rad = r( 9n ^- 3 \A 4 sin ( A 5 ^ (3 " +1)/2 ) + A 6 cos (A b t-^ n+1 ^ 2 ) 



(17) 
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(a) (b) 

FIG. 1: (a) R = -Rgr + SR as predicted in ref. [3l|] for the HSS model. The solid lines are the upper and lower envelopes of 
the solution, and the dashed line is -Rgr. We note the turning point in the lower envelope, at which point the amplitude of SR 
can become larger than -Rgr; (b) R = Rgr + SR for the AB model, where once again we observe a turning point in the lower 
envelope. 



where A\, A2, A3, A4, A$,Aq are constants. 

We have solved equation (fTS")) numerically, and R = Rgr + SR is shown in fig 1(a) taking a pure matter 
era as an example, so -Rgr = 4/3t 2 . We have used time coordinates such that -R vac = 1, and have chosen 
e = ((c2/2ci) 2ra /2c2) 1 ^ 2n+1 ^ -R vac = 0.1 and n = 1. We have evolved backwards in the time coordinate, over the 
range t = (0.25,0.1), using the initial conditions SR(ti) — 0.1, SR(ti) — 0, where ij = 0.25. By solving equation (fl5|) 
numerically, we have found that SR oscillates symmetrically around zero, confirming the behaviour expected from the 
WKB solutions. However, in fig 1(a) we have not exhibited the actual oscillations of R explicitly, since they are of too 
high freque ncy t o be resolved, and so have only presented Rgr and the upper and lower envelopes of the oscillations. 

From fig 1(a) we see that the lower envelope appears to have a turning point, indicating that in this approach 
after this time the oscillatory component SR will come to dominate over Rgr, and hence the Ricci scalar will become 
negative at some point in the past. 



B. AB model 



In refs. [32|, 123], a similar analysis has been considered for the model F(R) = R — i? vac /2 + xab- If we use xab in 
(fT4")) . we find the following equation 



S R + Uh-™2*\ SR + ^SR*0, 



which has the WKB solution 



SR* 



during the matter era and 



exp(l/ei 2 



Si sin (^J exp(2/3ei 2 )d^ + B 2 cos exp(2/3et 2 )dt 



(18) 



(19) 



SR ra a = ^ 3/4 exp(3a/4ei 3/2 ) 



B 3 sin ^ exp(a/2et 3/2 )d?j + B 4 cos (J exp(a/2et 3/2 )dt 



(20) 



for the radiation era, where a, B\, B2, -B3, -B4 are constants. We have solved equation (|18[) numerically, taking -Rgr - 
4/3t 2 , using time coordinates with R vac = 1 and choosing e = 0.32. The resulting -R = Rgr + SR is shown in fig 1(b) 



using the initial conditions SR(0A2) = 0.01, SR(0A2) — over the time range t = (0.42, 0.375). Once again, we have 
shown only the envelopes of the oscillations of R. These solutions exhibit similar oscillatory and growing behaviour 
as found for the HSS model. 
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C. Issues with the perturbative analysis 

Using this approach for both the HSS and AB models, we find similar behaviour for R. Specifically, SR undergoes 
rapid oscillations, and both the amplitude and frequency of these oscillations increases to the past. Since the amplitude 
of 5R grows without bound as t — > 0, it follows that SR will violate the condition SR <C -Rgr a t some point in the 
past, at which point the perturbative analysis will break down. Beyond this, we cannot assume that R = i?GR + SR 
is a solution to the gravitational field equations. 

This oscillatory behaviour presents us with a number of problems. The first is that since SR grows to the past, it 
will eventually satisfy \SR\ ~ Rgr- Beyond this point, R will periodically be negative, and when R is negative we can 
no longer use the expansion F(R) ss R — R vac /2 + x(-R). This is a problem particularly relevant for the HSS model, 
since for R < there will be a point at which F"(R) — 0, which is a singular point in the field equations. The second 
problem associated with the oscillatory behaviour of SR is that the frequency grows to the past without bound. It 
has been pointed out in ref.[3l| that this issue can be ameliorated by introducing an additional term oc R 2 /M 2 into 
the action, where M is a mass scale. This will provide a cut off in the frequency growth at w = M. 




i t 
(c) (d) 

FIG. 2: (a) the envelope functions for R in the HSS model, with e = 0.1, obtained by solving the full field equations numerically 
as described in the text, with perturbed initial conditions. We note that R differs significantly from the linearized approximation 
R — Rgr. + SR, obtained in the previous section; (b) SH = (H — Hgr)/H. It is clear that H oscillates, and the amplitude of 
these oscillations grows to the past; (c) the deviation of the scale factor from its General Relativistic limit, 8a = (a — a,GR,)/a. 
We see that a will deviate from agr as we evolve backwards in time, however Sa is highly suppressed compared to the deviations 
in R; (d) Sa over a smaller range of t close to the end point of the evolution. This shows the oscillatory behaviour of the scale 
factor which eventually develops. 
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(c) (d) 

FIG. 3: (a) the envelopes of SR = (R — Rgr)/R for the HSS model, obtained by solving the full field equations numerically, 
with unperturbed initial conditions. Since 7? = Rgr is not a solution to the field equations, we find that SR 0; (b) 
SH — (H — Hgk)/H. We see that H oscillates, and the amplitude of these oscillations grows to the past. Further, SH does 
not oscillate around zero, indicating that the Hubble parameter will deviate from Fgr as we evolve backwards in time; (c) SH 
over a small time regime to explicitly show the oscillatory behaviour of H; (d) 8a — (a — aGit)/ffl. We see that a will deviate 
from ogr as we evolve backwards in time, however 8a is highly suppressed. 



IV. NON-PERTURBATIVE NUMERICAL SOLUTION FOR THE RICCI SCALAR 

In this section, we first show that the above linearized analysis is only valid over a very limited range of R, and 
that the gravitational field equations are generically non-linear. To see this, we first write (0 as 

. y>" . „ R+ T/M$ Ry' - 2y 
R+SHR+^R 2 + pl - U \ „ X W 0, (21) 

where we have expanded the function F(R) w R — R vac /2 + x(R), an d neglected the term i? vac /2. Both of these 
assumptions are acceptable if we consider a regime where R > i? V ac- In section Hill we linearized this equation in SR, 
and in doing so we neglected the term (F m / F")(5R) 2 , which was assumed to be second order. However, since SR 
oscillates with high frequency w, it follows that the term (SR) 2 will not remain small, since it will grow like oj 2 to the 
past. Hence the accuracy of the above linearized analysis will begin to reduce when (F" 1 / F")(SR) 2 is of the same 
order as 3HSR. Beyond this point, we can no longer use the linearized equation for SR, and we must include these 
non-linear terms. 

For the AB model, we can explicitly show the range of SR over which the linearized analysis is valid. To see this, 
we consider the expansions ([SlfTU]). As an example, we take F'(R), which was expanded as F'(R) fts F'(Rcr) + 
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t t 
(c) (d) 

FIG. 4: (a) the envelope functions for R in the AB model, with e = 0.32, obtained by solving the full field equations numerically 
as described in the text, with perturbed initial conditions. As in the HSS model, we see that R differs significantly from the 
linearized approximation R = Rgh + SR, obtained in the previous section; (b) SH = (H — Hgr)/H, which oscillates (not 
around SH = 0); (c) Sa — (a — hgr)/b. We see that a will deviate from ogh as we evolve backwards in time, however Sa is 
highly suppressed; (d) Sa over a smaller range of t, which explicitly shows the oscillatory behaviour of the scale factor. 



F"(R GR )SR. For the AB model this reads F'(R) m 1 - e ~ flGR / e + \e~ RGR / e 8R. We compare this to our actual 
function F 1 (Rgr + SR) , given by F'(Rg R + SR) w i_ e --Rc!RA e -'5#A ) a nd we see that the expansion of F 1 (Rqk + SR) 
is only valid in the very limited range SR <C e — R vac /^b. Once this condition is violated, we can no longer use the 
linearized equation (jTTJ) for SR. Since in general we will have -Rgr 3> e — ^?vac/4fo, then it follows that the linearized 
analysis will break down for the AB model long before SR ~ Rgr- 

The above reasoning suggests that although the linearized analysis successfully predicts the oscillatory behaviour 
of these models, we can only trust the solution for a limited range of R, and beyond this the field equation for R is 
inherently non-linear. Wc conclude that we cannot find a solution to (JSJ) for which R ss Rgr f° r a ll R- Therefore in 
this section, we do not consider solutions to the field equations that are perturbations around the General Relativistic 
Ricci scalar, but rather look for full numerical solutions to |5]). However, it is not just a second order differential 
equation for i?, but rather a fourth order non- linear equation for the scale factor a(t). Therefore we will treat l[5]). 
H + 2H 2 = R/6 and H = a/a as a set of coupled non-linear differential equations for R, H and a. We stress that 
we are not making the assumption that R ~ Rgr + SR, and the only assumption that we will make is that during 
the matter era Rgr 3> Rvac- We do so because we wish to compare the full numerical results obtained here with 
the results of the linearized analysis of the section IIIII where the vacuum curvature was neglected. We have checked 
that introducing the vacuum curvature does not significantly effect our results, since we can absorb it into the Energy 
momentum tensor. 

We look for solutions to these equations for two sets of initial conditions. First we perturb R initially from its 
General Relativistic limit, and compare the full numerical solution to the linearized analysis of the previous section. 
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FIG. 5: (a) envelope functions for 5R = (R — Rgr)/R for the AB model, with unperturbed initial conditions. Since R = Rgb. 
is not a solution to the field equations, we find that 5R oscillates; (b) SH = (H — HgyC)/H. It is clear that H oscillates, and 
the amplitude of these oscillations grows to the past. Further, SH does not oscillate around zero, indicating that the Hubble 
parameter will deviate from Hgh as we evolve backwards in time; (c) SH over a small time regime, explicitly showing the 
oscillatory behaviour of H; (d) Sa — (a — aGR,)/a. We see that a will deviate from agr as we evolve backwards in time, however 
Sa is also highly suppressed. 

We then set the initial conditions such that R, a and H are initially exactly their General Relativistic values, and 
observe how the HSS and AB models evolve backwards through the matter era. 

A. HSS model 

For the HSS model, using the example of a pure matter era for which Rgr = —T/M^ = 4/3i 2 , wc find 



R + 3HR-2(n 



' R 



6n(2n+ l)e 



2n+l 



R 



■ T/M? 



pi. 



R 



2n+2 



(n + 1) 
3n(2n + 1 



-R 2 = 0, 



(22) 



together with the equations for H and a. The trace of the energy-momentum tensor T is related to a and H through 



the conservation equation V a T au = 0. Specifically, during the matter era we have T/M 



pi 



oc a 



To derive J22J from 



(0), we have used xhss in (|2Tj) . The only assumptions made are that F(R) w R — i? vac /2 + Xhss and that ignoring 
the i? V ac term will not significantly effect our results, which is valid throughout our numerical calculation. It is clear 
that R = —T/M^ is not a solution of these equations, and hence the HSS model will not exactly mimic General 
Relativity. As we will see, even if we set the initial conditions such that i?, H and a are initially equal to their General 



Relativistic values, if we evolve this model backwards in time then they will deviate from Rqr = —T/M 



P l- 
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We have solved the coupled differential equations ([22]) . H + 2H 2 — R/6 and H — a/a numerically over the range 
t = (0.25, 0.1), taking n = 1, R vac = 1 and e = 0.1 as before. We have found that we can only evolve the Ricci scalar 
over a small range of t. As we will see in the next section, this is because the Ricci scalar will generically evolve to a 
singularity at some finite point in the past. For now, we will ignore this singular behaviour and solve (|22[) for R over a 
small dynamical range. The Ricci scalar, obtained using a stiff differential equation solver, is shown in fig 2(a)[ Here 



we have chosen the initial conditions R{U) = (4/3t?) + 0.1, R(U) = -8/3if, H(U) = 2/3U, and a(ti) = (3/4) 1 / 3 tf /3 , 
where t\ = 0.25, and evolved backwards in the time coordinate. These initial conditions are the same as those chosen 
in the previous section, so that we can compare the full solution found here to the solution found using the linearized 
approximation R = Rgr + SR. In fig 2(a) we have exhibited the envelope of the oscillations of the Ricci scalar, 
and the General Relativistic solution Rqr- From the envelope functions, we see asymmetric oscillations of R about 
Rgr- We note that the solution obtained in this section is significantly different to the one obtained in the linearized 
analysis, shown in fig 1(a) Specifically, in fig 2(a) we find that there is no turning point in the lower envelope in the 
time range considered. Conversely, we find that the upper envelope function increases at a faster rate than predicted 
in the linearized analysis. 

We have also plotted SH = (H — Hqr)/H for this model in fig |2(b)| which is the fractional difference between 
the Hubble parameter H for the HSS model and the General Relativistic Hubble parameter during the matter era, 
Hgr = 2/3i. We see that 8H oscillates asymmetrically and not exactly around its General Relativistic limit 5H = 0. 
The amplitude of these asymmetric oscillations is highly suppressed but increasing to the past, and unless e is chosen 
to be sufficiently small these oscillations may come to dominate H. For the scale factor, we have again plotted the 
fractional difference 8a = (a — aQ^)/a in figs 2(c) and 2(d) and we see that a deviates from the General Relativistic 
scale factor, but this deviation 5a is highly suppressed. However, as with the Hubble parameter, 5a grows to the past, 
and hence may become significant at some earlier time. 

Having solved the system of equations for R, H and a by perturbing R initially from the General Relativistic 
limit, we now solve with no initial perturbation, that is we take the initial conditions R(U) = 4/3tp, R(t{) = -8/3tf, 
H(ti) = 2/3tj, and a(ti) — (3/4) 1 / 3 i?^ 3 , where ti = 0.25. The solution to equation (f2"2")) is shown in figs.3(a-d). Again 
we have taken e = 0.1, R vac = 1 and evolved over the time regime t = (0.25,0.1). We see that R oscillates, and H 
and a deviate from their General Relativistic values. However, once again 5a and 5H are suppressed, and this model 
closely mimics General Relativity over this range of time. 



B. AB model 

We can perform similar calculations for the AB model. Taking once again Rqr — 4/3i 2 as a specific example, we 
find 

R + 3HR - ?- + | (R + T/M 2 ) e R / e + ^-(^+2) =0. (23) 



We obtain (|23|) by using xab in equation (|2Tj) . We have solved this equation numerically, along with the equations 
for H and a, using the same differential equation solver as for the HSS model, taking e = 0.32 and using the initial 

conditions R(ti) = 4/3i 2 + 0.01, R(U) = -8/3if, H{U) = 2/3U and a{U) = (S/^/H^ 3 . As in sectionHlJ for the AB 
model we only solve equation (|23p over the very small time range t = (0.42, 0.375), due to the presence of a singularity. 
By solving (|23p over this limited dynamical range, we find asymmetric oscillations of R, H and a. The envelope of the 



oscillations of R are shown in fig 4(a) Once again, we see no turning point in the lower envelope, suggesting that R 
does not become negative for these models. We also observe that the upper envelope grows faster than predicted by 
the linear analysis. The Hubble parameter and scale factor undergo asymmetric, suppressed oscillations which grow 
to the past, as seen in figs.4(b-d). 

We also solve equation (|2"3")) without perturbing R, H and a from their General Relativistic limits initially, and 
evolving the system backwards over the same time regime t = (0.42, 0.375). The results are shown in figs.5(a-d). We 
find that R oscillates, and 5H = (H — Hqr)/H and 5a = (a — acaj/a grow, indicating that a and H both diverge 
from their General Relativistic limits. However, this divergence is highly suppressed, as was found in the HSS model. 



V. IMPROVED PERTURB ATIVE APPROACH 



So far, we have considered two approaches to solving the modified gravitational field equations. It has been found 
that the linearized approach, in which the ansatz R = Rqr + SR with 5R <C Rgr is used in (0, will not give a 
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solution that is indicative of the full solution, since non-linear terms cannot be neglected due to the rapid oscillations 
of the Ricci scalar. The second approach considered was to solve the full gravitational field equations numerically. In 
doing so, it was found that the Ricci scalar undergoes non-linear oscillations, very closely (but not exactly) about its 
General Relativistic limit. However, we found that solving the gravitational field equations over significant timescales 
is impossible due to some kind of singularity at a finite time. 

In this section, we consider neither a linearized analysis nor a full numerical study of the field equations. Instead, 
under sensible assumptions, we show that the trace of the gravitational field equations for both the AB and HSS 
models reduces to a non-linear wave equation, and hence that R undergoes asymmetric oscillations about its General 
Relativistic limit, in agreement with the results of section HVl Moreover we find that there is a singularity in R at a 
finite time in both models. Since the approach taken is model specific, we consider the AB and HSS models separately. 

A. AB model 

We begin by substituting %ab into ([5]), giving 

- iUe~ R / e -R + i? vac -(R+ 2e)e- R l f - = — (24) 



Next, we define z = e R ^ e , so R — — elogz and (|24p becomes 



T 



30s + elogz + i? vac +ez (log z- 2) = — j., (25) 



where we have assumed that i? vac "C TjM\ and hence neglected the vacuum energy. Next, by defining Z\ = ze T / cM pi, 



we can write (1231) as 



3Dzie T / eM p> + e log z x + ee T ' eM ^z x ( log z x + — ^ - 2 J =0. (26) 



To solve this equation, we write it as 



3e T/2eM^ ± ( £ T/2eM^ dfl Jl_ T/2eM^ dT\ , 3 ( T/2eM^ dfl , _£l T/2 e Af£ *[\ f^^dT 

6e dt[ 6 dt + eM^ B dt) +6 [ e dt + eM^ 6 dt ) 2eM p 2 , dt + [Z7) 

9e T/2eM^ a -lda f T/KM^dZl J* e T/2eM^\ + e]fJ + ^T/eM^ L + J_ _ ^ = q 

dt \ dt eM p 2 , dt J 6 \ eM^ J 

If we now introduce the fast time coordinate 

X=f e- T / 2M *dt, (28) 
which is the timescale associated with the oscillations of the Ricci scalar, then we can write (|2"Tf as 



(29) 

where primes denote derivatives with respect to A, and we have defined H — a' /a. If we solve equation (|29[) for z%, 
then we can deduce R from the relation 
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R 



T 

Mix 



e log z\ . 



(30) 



Note that (f3"0"|) is of the form R = -Rgr + -Rose, however we have not assumed that R osc R&R. to deduce this 
expression. 

If we now define x = Z\ — 1, then we can write (|29[) as 



3a;" 



2eM* x 



T' + 9H] x' + e log(l + x) + x 



2e*M* 



T 



rpfi 



HT' + ee T /^ 



- 2 



(31) 



3 T/eM 



p'(l +ar)log(l + .t) = - 



2e*M* 



T 



HT' 



T/cMp[ 



T 



- 2 



Equation (|3 1 [) is a non-linear, second order inhomogeneous differential equation for x, and the full solution is the sum 
of the particular solution to (f3"Tj) and the solution to the corresponding homogeneous equation given by setting the 
right hand side to zero, x = is now the General Relativistic limit, but it is clear that this is not an exact solution. 
We will find that x undergoes asymmetric oscillations governed by the homogeneous equation, around the particular 
solution of (l3Tj) . We now calculate the oscillatory and drift components of x. 



1. Drift 

We begin by calculating the drift of x away from x — 0. To do so, we make the assumption that the derivative terms 
x" and x' in (|3"Tj) are subdominant and hence can be neglected. With this assumption, (|3"Tj) reduces to an algebraic 
equation for x, 



elog(l + x) + x 



2e 2 M p 4 



rpfl 



HT' + ee 



T/eM* 



T 



+ ee T / eM pi (1 + x) log(l + x) (32) 



T 



rpfi 



HT' 



T 



- 2 



If we assume that then by expanding log(l + x) in (|32p we find the following expression for x, 



= - 3e T / £M p 



64 
9^ 



T 



3e 2 t 4 3 \eM^ 



- 2 



(33) 



This solution is valid if we only consider terms linear in x, that is we consider terms of order e T ^ EM p' <C 1 only. Taking 
the derivative of x, we find that x' and x" are given by 



e T/2eM^ 1± 



It T/eM 

x = e ' p 



T 



x + x 



(34) 
(35) 



Since x' and x" are much smaller than x (they are suppressed by a factor of exp(T/2eMp[) <C 1 and exp(T / eM^) <C 1 
relative to x respectively), then it follows that our original assumption that we can neglect derivative terms of x in 
(f3"Tj) is valid, and (f3"3"| is an approximate solution to (|31[) . 
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2. Oscillations 



The solution to the homogeneous equation, obtained by setting the left hand side of (f3Tj) to zero, will give the 
oscillatory component of x and hence R. We will first show that the last two terms on the left hand side of (|3"3")) , are 
highly suppressed, and hence can be neglected for T/AA e. To do so, we write T", (T') 2 and HT' in terms of t, 



IT = ( ? \ 2 = _j_ eT/eMa {9et 2 4) ^ = i6 r/e< 

JVfpj U 3 e- T/2eM pV ' AfJ 9i 6 e 1 h M p 2 , 9i 4 ' 1 j 



T 

and the corresponding homogeneous equation to ((33)) can be approximately written as 



3x"+ — —jT' + 9H x'+elog{l+x)+3e T/eM ^x 



2eM p 2 , 



64 8 1 



9e 2 t e Set 4 3 



+ -e log(l + x) + 



T 
^5 



+ee T/6M| 1 l og (l +a; ) = 0. 

(37) 



We see that the last two terms on the left hand side of (I37| have a common factor of exp(T / eM^) <C 1 and hence 
can be neglected. We arrive at the following equation for x, 



x" + f3x' + \e\og{l + x) w0, (38) 

where = „ |p T / + 3-ff is a small function of A. In ([38]) , we will assume that iJ = Hgr and T = Tqr, since 

the oscillations of H and a are suppressed by a factor of exp(T / eM^) and can therefore be neglected. With this 
assumption, (|38p reduces to an equation purely in terms of x and its derivatives. 

Equation (|38p describes an anharmonic, non-linear oscillator with a small negative damping term, and hence we 
have shown analytically that the Ricci scalar will undergo asymmetric oscillations, in agreement with the results of 
section HVl We henceforth describe the oscillator considered here as the logarithmic oscillator. The potential V{x) 
for the logarithmic and simple harmonic oscillators are shown in figs.6(a,c). The logarithmic potential will give rise 
to wave solutions that are not symmetric. 

Before continuing, we observe that if we expand R as R = Rqr + 5R for SR -C e, then x can be expanded as 
x = cxp(— (R — i?GR.)A) — 1 ~ —SR/e, and ([38]) becomes 



8R" + 138 R' + -e8R^0. (39) 

If we write (|39)) in terms of t, then it reduces to the linearized equation (fT8| . 

We have solved the non-linear equation ([38]) numerically, using the initial conditions x\ = exp(— 0.01/0.32), x[ = 0, 
and evolved over the range A = (—1350, —13339). This range of A and choice of initial conditions are the same as 
were chosen in section HVl where the full gravitational field equations were solved numerically. In figiZUa), we compare 
the Ricci scalar R — i?GR — elog(l + x) obtained in this section to that found by solving the full gravitational field 
equations numerically in section HVl It is clear that the results obtained in this section closely mimic the full numerical 
solution. In figlUb) we compare SR a = R—Rqr = — 6 log(l + x osc +Xdr) to 5R = R— Rgr, where R is the Ricci scalar 
obtained by solving the full gravitational field equations numerically. We see that the fractional difference between 
5R and 5R a remains small throughout the dynamical range considered. 



3. Existence of Singularity 



Using this approach, we have found that we can only evolve x and hence R backwards in the time coordinate t over 
the s ame l imited dynamical range as in section Hvl However, by solving equation (|38|) we can see why this is the case. 
In fig 7(c) we have perturbed x from x = and evolved backwards in the A coordinate. We see that the amplitude of 
the oscillations of x increase due to the presence of the damping term in (|38|) , and after a finite time x — > — 1. Since 
R = — — e log(l + x) and R — T/M^ — ex/(l + x), then it follows that both R and R will be singular at this point. 
We conclude that wide ranges of initial conditions x\ and x\ will give rise to a singularity in the field equations as we 
evolve backwards in time. That is not to say that all initial conditions will lead to this problem, and there exists an 
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extremely restricted range of x\ and x\ at the beginning of the matter era that will not give rise to singular behaviour 
at some finite time in the past. 

Instead of evolving backwards over the matter era, if we choose cc; appropriately at the beginning of the matter era 
and evolve forwards in time, x (and hence R) will be well behaved and regular. To find the allowed initial conditions, 
we can first use the fact that x\ > ~ 1 is the lower bound. To obtain an upper bound, we assume that the damping 
term in f38|) is negligible and consider 

.T" + ielog(l + a;)«0. (40) 
This expression can be multiplied by x' and integrated to obtain the following expression, 

\{x') 2 + ie(l + x) [log(l +x)-l] = E, (41) 

where E is a constant of integration. The values of x at the peaks and troughs of its wavetrain can be obtained from 
(|4"T) by setting x' — and solving the algebraic equation 



E = -e(l + x) [log(l +x)-l]. (42) 

For a given E, there will be two solutions for x, corresponding to the maxima and minima of the wave. We have 
deduced that x = — 1 is the minimum value that x can take in order to avoid singularities, and substituting this into 
(l42)l gives E = 0. As we have stated, (|42j) will yield two solutions for any particular E, and so finding the second 
solution for E = will give us the maximum value of x. This solution is x = exp(l) — 1, and hence we conclude that 
in order to obtain a matter era free from singularities, we must impose the condition that —1 < x < exp(l) — 1 at the 
beginning of the matter era. In terms of SR = R — Ron, this corresponds to the range — e < 5R < oo. 

Since we will generically encounter a singularity if we evolve x backwards in the time coordinate t, we now solve 
equation (|38|) by taking x initially at some early time and evolving forwards in the A coordinate, which corresponds to 
evolving forwards in the standard time coordinate t. In doing so, x and hence R will remain regular. We take as initial 
conditions x\ — exp(— 2.6/0.32) — 1, x[ = and evolve over the range A = (—13339,220). By solving (|3"5|) , we obtain 
x osc , which is the oscillatory component of x. The full solution to (j3Tj) is approximately given by X ' •Z'dr I *^osc j 

where 

Xdr is the drift term given in (|3"3"|) . In fig[5](a-d) we compare 5R a to 5R. We see that there are two competing effects; 
the non-linear oscillations, which are initially large but quickly decay, and the drift away from SR = 0, which grows 
as we evolve forwards in the time coordinate. The approximate solution obtained in this section closely mimics the 
full numerical solution. However, the predicted drift begins to deviate from the full solution at late times. However, 
this only occurs when the approximation T/M^ 3> e is no longer valid. 

4- Constant Energy Momentum Tensor 

As an aside, we conclude this section by considering a constant energy momentum tensor, that is we take Rqr = 
R a = const ;» e. For this Ricci scalar, we have A = exp(i? /2e)t, T' = and H = ^/i? /12 cxp(— R a /2e). Using these 
in (|29[) . we obtain the following equation for z%, 

z'l + 3 v / R /12e~ R °/ 2e z' 1 + \e log z x « 0. (43) 

The case of constant Rqr was discussed in ref. [33| using the perturbative analysis of section IIII| and it was found 
that the Ricci scalar contained an exponentially growing component, suggesting an instability in this model. Here, we 
solve (|43p for z\ and obtain the Ricci scalar from R — Rq — elogzi. The solution is presented in figsHJa-b), where we 
have taken e = 0.4, i? V ac = 1, Rq = 14i? vac and evolved over the time regime A = (0, 500). We see that R undergoes 
damped, rapid oscillations about its General Relativistic limit. We find no exponential growth, contrary to the claim 
made in ref. [33j . 
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FIG. 6: (a) the logarithmic oscillator potential Vi = — x + xlogx; (b) the oscillator potential in the HSS model, V2 = 
-(1 + z) 2 / 3 (3/2- (l + z) 1 / 3 ); (c) the simple harmonic oscillator potential. The asymmetry of the potential in (a) and (b) 
will lead to asymmetric solutions. 



B. HSS Model 



In this section, we consider a similar calculation to that presented above, but now for the HSS model. We first 
re- write the trace of the gravitational field equations in terms of a dimensionless function and a fast time coordinate. 
We then derive approximate expressions to describe the oscillations and drift of the Ricci scalar away from its General 
Relativistic limit. We begin by substituting xhss into (0, which gives 



M _2(n+l) (vi?)2 
R 



R 



T 



R 



2n+2 



n + 1 



M,2 J 6n(2n + l)e 2rl + 1 3n(2n + 1) 



Bf w 0. 



(44) 



where we have neglected the term i? vac . To transform equation (|44[) into an oscillator equation, we must make a series 
of field redefinitions. By first defining 



(-T/M2 



2 \2n+l 



R 2n+1 



>0, 



a(t) 



,2n+l 



i-T/M^) 



2 \2n+l 



>o, 



(45) 



we write (14411 as 



i 



,l/(2n+l) 



1 T 



in Ml 



au 



2n/(2rt+l) 



(46) 
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9 7 




(c) 

FIG. 7: (a) The Ricci scalar obtained by solving the full gravitational field equations for the AB model numerically (solid line) 
is exhibited along with the approximate solution R — Rgr. — elog(l + x) (dashed line), where a; is a solution to (|38[) . We see 
that the approximate solution closely mimics the full numerical solution; (b) The difference (SR a — SR)/SR, as defined in the 
text. We see close agreement between the full and approximate solutions; (c) x for the AB model (that is, the solution to (1381 1). 
We see that after a finite time x — > — 1, at which time R — > oo. This singularity generically occurs when we evolve the AB 
model backwards through the matter era. 

Next, by introducing a fast time coordinate 

(\ n+l/2 
if J * m 

and denning x = u — 1, we can write (|46[) in terms of u and A as 

*" + (igSW + 3 ") •> - s - + '^k a(l + If """' +1 ' (48 ' 

where primes denote differentiation with respect to A, and H = a' /a. This is a non-linear, inhomogeneous second 
order differential equation for x, and by solving for x we obtain the Ricci scalar from the relation R = —T/(M^{1 + 

x) 1 /( 2n+1 )). As before, we can deduce both the oscillatory component of R and the drift away from i?GR- 
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(d) 



FIG. 8: (a) We take x at an early time initially and evolve forwards in the time coordinate A (and hence forwards in t). The 
solid dark line is SR = R — Rgr, obtained by solving the full field equations numerically, and the light dashed line is the 
approximate solution <5i? a = — elog(a;osc +a;dr), where a; osc is the oscillatory component of x, given by the solution to (|38|l . and 
Xdr is given in (|33|l : (b) The full and approximate solutions SR and SR a over a small time regime. We see close agreement 
between the two; (c) SR and SR a over a different time regime, showing the drift away from SR = 0. Again we see a close 
agreement, except at late times where R ~ i? V a C ; (d) The difference 8R a — SR. This difference is small for Rgr 3> e, but 
increases when -Rgr ~ Rvmc- 





04 0.6 0E 



1 1.2 1.4 IB 1.1 



(a) 



(b) 



FIG. 9: (a) R — Ro — elogzi, for constant _Ro = 14. We see that the Ricci scalar undergoes damped oscillations around Rq; (b) 
is a phase diagram of yi = ii against z\. We see that zi undergoes damped oscillations around zi = 1, which is the General 
Relativistic limit. 
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FIG. 10: (a) The Ricci scalar obtained by solving the full gravitational field equations for the HSS model numerically (solid 
line) is exhibited along with the approximate solution R = Rgr/ (1 + x) 1 ^ 3 (dashed line), where a: is a solution to (|5ip . We see 
that the approximate solution closely mimics the full numerical solution; (b) The difference (x a — x)/x, as defined in the text. 
We see close agreement between the full and approximate solutions; (c) x for the AB model (that is, the solution to (J38J ) . We 
see that after a finite time x — > —1, at which time R — ► oo. This singularity generically occurs when we evolve the HSS model 
backwards through the matter era. 

We now take as a specific example a pure matter era and set n = 1, in which case we can write a — £ 3 /-Rgr and 
A = — 16/27e 2 i 3 . We also make the assumption that H rss Hgr, which implies that 

/ '-^'-K(TT^- 1 ) + 3^ (1 + a:> " ' (49) 

where we have neglected the third term on the left hand side of (|48|) . which is negligible. 



1. Drift 



As in the AB model, an approximate expression for the drift can be found by solving the full inhomogeneous 
equation (|49]) . We take as an ansatz 



with this choice, x' and x" are of order x' ~ A -3 and x" ~ A -4 , and hence the derivative terms in P5|) can be 
neglected if we only consider terms of order O (A -2 ). As a result, it reduces to an algebraic expression for x, which 
is solved by (this can be verified by direct substitution). 
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(c) (d) 

FIG. 11: (a) We take x at an early time initially and evolve forwards in the time coordinate A (and hence forwards in t). 
The solid dark line is x = (Rgr/R) 3 , obtained by solving the full field equations numerically, and the light dashed line is the 
approximate solution x a = Kosc + Xdr, where x osc is the oscillatory component of x, given by the solution to (|51[) . and Xdr is 
given in (|50[) ; (b) The full and approximate solutions x and x a over a small time regime. The oscillatory components of x and 
x a show close agreement oscillations; (c) x and x a over a different dynamical range. We see the drift away from x — 1; (d) The 
difference x — x a . This difference is small for -Rgr 3> e, but diverges when -Rgr ~ Rvac- 

2. Oscillations 

Having calculated the drift by solving the oscillatory component of x can be obtained by neglecting the last 
term on the right hand side of (|49p (since it is of order ~ A~ 2 ), and solving the equation 

» 10 , e / 1 \ 

* -^-HOT^- 1 )" - (51) 

This equation describes the non-linear oscillations of x around x — 0, that is around the General Relativistic limit 
R = — T/Mpj. We have solved equation (f5Tj) . using the initial conditions X\ = (-Rgr/(^gr + 0-1)) 3 ~ 1) x i — 0j over 
the range A = (—3790,-69260) and taking e = 0.1, and the results are exhibited in fig. (fTT)]) . We compare the Ricci 
scalar obtained in this section, given by R — i?GR./(l + a^) 1 ^ 3 , to the Ricci scalar found in section UVl and find that 
they are in close agreement. 

3. Existence of Singularity 



As in the AB model, by using this method we find that we can only evolve x backwards over a very limited 
dynamical range, and this is once again due to the presence of a singularity. In fig 10(c) we see that as we evolve x, 
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after a finite time x — > —1, and since i? = -Rgr/(1 + x) 1 / 3 , we find that as x — > — 1, i? — ► oo. This singularity arises 
due to the initial conditions that we have imposed. 

As in the AB model, to evade this singularity we will instead choose initial conditions at an early time and evolve 
forwards through the matter era. By doing so, x will remain regular throughout. To find the initial conditions that 
x may take in order to evade this singularity, we perform the same steps as in the AB model. We first observe that 
x; n > — 1 is a lower bound, and the upper bound can be obtained from the expression 



M! _ (1 + x yn/(ln+l) ( ^±1 _r 1 + x) l/(2n+l)\ = E (52) 

2 \ 2n J 

where we have neglected the damping term and terms of order ~ A -2 in (|5~lj) . By substituting x — > —1, we find 
that E — 0, and hence the upper bound for x\ is found by solving ()52|) for E = 0. Doing so, we arrive at Xi < 
((2tt, + l)/2n) 2n+1 - 1, which for n = 1 is given by x; = (3/2) 3 - 1 = 2.375. In terms of R = Rqr/(1 + x) 1/3 , this 
corresponds to the choice 2nRQ^/(2n + 1) < i?; n < oo at the beginning of the matter era. 

With this in mind, we now solve (|5ip by evolving forwards in the fast time coordinate A, which corresponds to 
evolving forwards in the t coordinate. We solve (|5ip over the range A = (—69260,-475), with initial conditions 
X[ = (^gr/(-^GR + 2.6)) 3 and x\ — 0. By solving (|5"Tj) we obtain the oscillatory component x osc of x, and hence 
x ~ x osc + Xdr, where x^ T is given in ([50]) . In fig. pdj) we compare x obtained in this section to x = (Rgr/R) 3 , where 
R is the Ricci scalar obtained solving the full field equations numerically. As in the AB model, there are two effects; 
the non-linear oscillations of x (and hence R) which decay as we evolve forwards in time, and a drift in x away from 
its General Relativistic limit x = 0. It is clear that the approximate solution obtained in this section closely mimics 
the full solution of section IIV1 



VI. DETERMINATION OF THE HUBBLE PARAMETER 

In section [V] we have calculated the Ricci scalar for the HSS and AB models, and have found that it will undergo 
non-linear oscillations and will drift away from the General Relativistic limit R = —T/M^ as one goes back in time. 

From the Ricci scalar we can now calculate the Hubble parameter for the two models, using H + 2H 2 — R/6. We 
look for a solution of the form H = i?GR + SH, where SH <C Hgr- A solution of this form must exist in order for 
these models be viable; they must have a matter era for which H w 2/3t for example, in order for normal structure 
formation to take place. Linearizing in 5H, we find the following expression 

±5H + 4H GR 6H= 5 -^. (53) 
Next, by using the time coordinate A, where A cx (i 7 "' ) -1 / 2 , we obtain 

where oqr is the General Relativistic scale factor. This expression can be integrated to obtain SH, 



SH =^ + -±- [ a M^-d\^^ + X fsRdX, (55) 
«gr gr J A a GR 6A J 

where C is a constant of integration, and we have taken the slowly varying factor Oq R /A outside the integral, which 
is a good approximation for A ^> 1. We see that SH contains two terms; one describing the oscillations of H due to 
8R, and a term that goes like SH ~ Ca^p. The oscillatory component is suppressed by a factor of A -1 , and hence 
the oscillations of R will not have a significant impact on the Hubble parameter. We note in obtaining SH, we have 
not had to specify either SR or A, and hence (|55|) is valid for any model for which RF"(R) <C 1, including the HSS 
and AB models. 

In a similar manner, we can also consider the scale factor. By writing H = a/a and expanding as a — acR + Sa, 
we have at linear order 
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u m « _ h GK . Sa a GR Sa 

H GK + 5H = -tt h 2 — • ( 56 ) 

a «GR OGR a GR 



Next, by introducing the fast time coordinate A, we write (jSij)) as 

d ( Sa \ 5H 



dX UgrJ " A ' (5?) 
and by using (]55f) and integrating, we arrive at the following approximate expression for 5a, 

f S II a f f 

5a w Da G R + a GR / -j-dA « Da G R + — 7- / <5MA + Ca G R / — ; — r , (58) 



A 6A 2 J J a GR A 



where D is a constant of integration. 5a possesses an oscillatory term due to the rapid oscillations of 6R, which is 
suppressed by a factor of A~ 2 <C 1. We conclude that the rapid oscillations of the Ricci scalar has no significant 
impact on the scale factor for models in which RF"(R) -C 1. 



VII. DISCUSSION AND CONCLUSIONS 



In this paper we have solved the gravitational field equations for the HSS and AB models numerically and by 
using an alternative perturbative analysis. We have shown that the oscillations of the Ricci scalar arc inherently 
non-linear, although they can be modeled as linear waves in a certain regime. However, for the AB model this linear 
regime is particularly restricted. Rather than using the earlier perturbative analysis, we have re-written the trace of 
the gravitational field equations as a damped, driven, non-linear oscillator. By solving this new equation, we have 
found that the Ricci scalar does indeed oscillate with high frequency, as predicted in ref.[3l|. As the amplitude and 
frequency of these waves grow to the past, we see non-linear behaviour becoming increasingly important. We have 
also found that the Ricci scalar does not exactly oscillate around its General Relativistic limit, but rather there is a 
highly suppressed drift. 

In section IPVl we solved the full gravitational field equations for both models numerically. By specifying the initial 
conditions at late times and evolving backwards, we observed that in general the Ricci scalar would evolve to a 
singularity after a finite time. This is another important consequence of the non-linear terms in the field equations, 
since no such behaviour was observed when using the linearized approach. By using our oscillator equation, we were 
able to explain why this singularity occurs, and found that it can be avoided by choosing initial conditions at some 
early time and evolving forwards. We were able to derive the allowed initial conditions for both the AB and HSS 
models. We have found that although the singularity could potentially be evaded with a suitable choice of initial 
conditions, the resulting Ricci scalar is likely to be unstable to perturbations away from a perfect matter era. To 
obtain a viable modified gravity model, it is clear that some method of regularizing this singularity is required. 

Finally, we have considered the effect of the oscillations of R on the Hubble parameter and scale factor. We have 
found that the oscillatory components of H and a are suppressed by factors of A -1 and A~ 2 respectively, where 
A cx \/F". We have concluded that 5H and 5a will remain small, in spite of the potentially large oscillations of R. 
This conclusion is generic to models which have RF" <C 1 • 

Although we studied specifically the HSS and AB models in this paper, we believe that many of the results obtained 
are generic to F(R) theories of gravity which satisfy RF" -C 1. We expect that the procedure adopted in section IVl 
that is writing the trace of the gravitational field equations as a non-linear oscillator equation, will generalize to all 
models for which RF" <C 1, and if that is the case then the singularity encountered in this paper will be a common 
feature of F(R) models. 
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